home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / fft.lha / fft / fft_input.cc < prev    next >
C/C++ Source or Header  |  1993-08-08  |  13KB  |  465 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *                  Fast Fourier Transform
  6.  *
  7.  *              Processing the Input Sequence
  8.  *
  9.  * The functions below handle different cases of the input (real/complex,
  10.  * with/without zero padding), form the array A, and perform the radix 2 FFT
  11.  * algorithm itself, taking advantage of the particular form of input.
  12.  *
  13.  * General radix 2 FFT algorithm is as follows
  14.  *    Input: x_re[i], x_im[i]        the input sequence
  15.  *    Output: Aj = SUM[ Xk * W^(k*j) ],  j,k = 0:N-1,
  16.  *    where
  17.  *        W = exp(-2pi I/N)
  18.  *        Xk = ( k<N/2 ? x_re[k] + I x_im[k] : 0.0 ) with zero padding
  19.  *        Xk = x_re[k] + I x_im[k], without zero padding
  20.  *                                
  21.  * Algorithm
  22.  *    1.  Fill in the complex array A performing the permutation
  23.  *        of input data
  24.  *        Ai = ( j < N/2 ? x_re[j] + I x_im[j] : 0.0 ), j = 0..N-1
  25.  *        or
  26.  *        Ai = ( x_re[j] + I x_im[j] ), j = 0..N-1
  27.  *        where 
  28.  *         index i is the "inverse" of the index j (in the sense that
  29.  *        m-bit string representing the value of 'i' is read from right
  30.  *        to left the string corresponding to index 'j'; m = log2(N)
  31.  *
  32.  *    2. Transform itself
  33.  *       n2_l = 1;            // 2^l at l=0
  34.  *       for l=0 to m-1 do        // m = log2(N)
  35.  *           for k=0 by 2* n2_l to N-1 do
  36.  *               for j=0 to n2_l-1 do
  37.  *            compl w = exp(- I j*pi/2^l );
  38.  *                i1 = j + k;
  39.  *            i2 = i1 + n2_l;
  40.  *            A[i1] = A[i1] + A[i2]*w        // Batterfly operation
  41.  *            A[i2] = A[i1] - A[i2]*w    
  42.  *           od
  43.  *             od
  44.  *           n2_l *= 2                // 2^l at the next l
  45.  *       od
  46.  *
  47.  * The first three steps of the algorithm above can be performed with
  48.  * simplified formulas from
  49.  *        G.Nussbaumer. Fast Fourier Transform and Convolution
  50.  *        Algorithms, M., Radio i sviaz, 1985, p.135 (in Russian)
  51.  * The formulas carry out the 8-point FFT for each set of 8 points
  52.  * A[k+0], A[k+1], ... A[k+7], k=0,8,16..N-8     
  53.  *
  54.  * The formulas are given below for references, x0-x7 being the source
  55.  * data, y0-y7 the transformed data. Moreover, one has to keep in mind
  56.  * that the source data x0-x7 have been already permuted (as opposite to
  57.  * the book that assumes the source data to be arranged in the natural
  58.  * order rather than permutted one)
  59.  *
  60.  *    t1 = x0 + x1       m0 = t7 + t8           s1 = m3 + m4    y0 = m0
  61.  *    t2 = x2 + x3       m1 = t7 - t8           s2 = m3 - m4    y1 = s1 + s3
  62.  *    t3 = x4 - x5       m2 = t1 - t2           s3 = m6 + m7    y2 = m2 + m5
  63.  *    t4 = x4 + x5       m3 = x0 - x1           s4 = m6 - m7    y3 = s2 - s4
  64.  *    t5 = x6 + x7       m4 = cu*(t3 - t6)                      y4 = m1
  65.  *    t6 = x6 - x7       m5 = I * (t5 - t4)                     y5 = s2 + s4
  66.  *    t7 = t1 + t2       m6 = I * (x3 - x2)                     y6 = m2 - m5
  67.  *    t8 = t4 + t5       m7 = -Isu* (t3 + t6)                   y7 = s1 - s3
  68.  *
  69.  *    cu = su = sin( pi/4 )
  70.  *
  71.  ************************************************************************
  72.  */
  73.  
  74. #include "fft.h"
  75. #pragma implementation "fft.h"
  76.  
  77. #include <math.h>
  78. #include <Complex.h>
  79. #include <std.h>
  80.  
  81. /*
  82.  *-----------------------------------------------------------------------
  83.  *            Perform the transform itself
  84.  * on the complex array A with inplace. First three steps of the algorithm
  85.  *            are assumed to have been already performed.
  86.  *
  87.  */
  88.  
  89. void FFT::complete_transform()
  90. {
  91.   register int n2_l = 8;        // 2^l at l=3
  92.   register int inc = N/2/n2_l;        // Increment for the index in W
  93.  
  94.   for(; n2_l < N; n2_l *= 2, inc /= 2)
  95.   {
  96.     register Complex * ak;
  97.     for(ak=A; ak < A_end; ak += 2*n2_l)
  98.     {
  99.       register Complex * aj1 = ak;
  100.       register Complex * aj2 = aj1 + n2_l;
  101.       register Complex * w = W + inc;    // Since j=0 case is handled explicitly
  102.  
  103.       Complex ai1  = *aj1;
  104.       Complex ai2w = *aj2;             // Batterfly operation in case of
  105.       *aj1++ += ai2w;            // j=0 ==> W=1
  106.       *aj2++ = ai1 - ai2w;
  107.       for(; aj1 < ak+n2_l; w += inc)
  108.       {
  109.         Complex ai1 = *aj1;        // w = exp(-I 2pi/N j*inc) =
  110.         Complex ai2w = *w * *aj2;       // exp(-I j pi/2^l)
  111.         *aj1++ += ai2w;
  112.         *aj2++ = ai1 - ai2w;        // Batterfly operation
  113.       }
  114.     }
  115.   }
  116. }
  117.  
  118. /*
  119.  *------------------------------------------------------------------------
  120.  *    The most general case of complex input without zero padding
  121.  *
  122.  * The Nussbaumer formulas are applied as they are, no further simplification
  123.  * seems possible.
  124.  */
  125.  
  126. void FFT::input(
  127.          const Vector& x_re,    // [0:N-1] vector - Re part of input
  128.          const Vector& x_im)    // [0:N-1] vector - Im part of input
  129. {
  130.   are_compatible(x_re,x_im);
  131.   if( x_re.q_lwb() != 0 || x_re.q_upb() != N-1 )
  132.     _error("Sorry, vector [%d:%d] cannot be processed.\n"
  133.        "Only [0:%d] vectors are valid",
  134.        x_re.q_lwb(), x_re.q_upb(), N-1);
  135.  
  136.   register int k;
  137.   register Complex * ak;
  138.  
  139.   for(ak=A,k=0; k < N; k +=8)
  140.   {
  141.     const double cu = W[ N/8 ].real();            // cos( pi/4 )
  142.  
  143.     register int i;
  144.     i = index_conversion_table[k];    // Index being "reverse" to k
  145.     Complex x0(x_re(i),x_im(i));
  146.  
  147.     i = index_conversion_table[k+1];
  148.     Complex x1(x_re(i),x_im(i));
  149.  
  150.     i = index_conversion_table[k+2];
  151.     Complex x2(x_re(i),x_im(i));
  152.  
  153.     i = index_conversion_table[k+3];
  154.     Complex x3(x_re(i),x_im(i));
  155.  
  156.     i = index_conversion_table[k+4];
  157.     Complex x4(x_re(i),x_im(i));
  158.  
  159.     i = index_conversion_table[k+5];
  160.     Complex x5(x_re(i),x_im(i));
  161.  
  162.     i = index_conversion_table[k+6];
  163.     Complex x6(x_re(i),x_im(i));
  164.  
  165.     i = index_conversion_table[k+7];
  166.     Complex x7(x_re(i),x_im(i));
  167.  
  168.     Complex t1 = x0 + x1;
  169.     Complex t2 = x2 + x3;
  170.     Complex t3 = x4 - x5;
  171.     Complex t4 = x4 + x5;
  172.     Complex t5 = x6 + x7;
  173.     Complex t6 = x6 - x7;
  174.     Complex t7 = t1 + t2;
  175.     Complex t8 = t4 + t5;
  176.  
  177. #define m2 t1
  178.     m2 -= t2;                // Forget t1 from now on
  179. #define m3 x0
  180.     m3 -= x1;                // Forget x0 from now on
  181.     Complex m7 = t3 + t6; m7 = Complex(cu*m7.imag(),-cu*m7.real());
  182. #define m4 t3                // Forget t3 from now on
  183.     m4 -= t6;
  184.     m4 *= cu;
  185.     Complex m5 = t5 - t4; m5 = Complex(-m5.imag(),m5.real());
  186.     Complex m6 = x3 - x2; m6 = Complex(-m6.imag(),m6.real());
  187.  
  188. #define s1 m3
  189.     Complex s2 = m3 - m4;        // Forget m3 from now on
  190.     s1 += m4;
  191. #define s3 m6                // Forget m6 from now on
  192.     Complex s4 = m6 - m7;
  193.     s3 += m7;
  194.  
  195.     *ak++ = t7 + t8;                 // y0  
  196.     *ak++ = s1 + s3;                // y1
  197.     *ak++ = m2 + m5;                // y2
  198.     *ak++ = s2 - s4;                // y3
  199.     *ak++ = t7 - t8;                // y4
  200.     *ak++ = s2 + s4;                // y5
  201.     *ak++ = m2 - m5;                // y6
  202.     *ak++ = s1 - s3;                // y7
  203.   }
  204.   assert( ak == A_end );
  205.   complete_transform();
  206. }
  207. #undef m2
  208. #undef m3
  209. #undef m4
  210. #undef s1
  211. #undef s3
  212.  
  213. /*
  214.  *------------------------------------------------------------------------
  215.  *           Real input sequence without zero padding
  216.  * When x[j]-s are all real, some intermediate results are either pure
  217.  * real, or pure imaginaire, which are much cheaper to compute than
  218.  * the complex ones. In Nussbauner's formulas above, t1 through t8,
  219.  * m0 through m4, s1 and s2 are all real, whilest m5 through m7, s3, and
  220.  * s4 are pure imaginaire.
  221.  *
  222.  */
  223.  
  224. void FFT::input(
  225.          const Vector& x_re)    // [0:N-1] vector - Re part of input
  226. {
  227.   x_re.is_valid();
  228.   if( x_re.q_lwb() != 0 || x_re.q_upb() != N-1 )
  229.     _error("Sorry, vector [%d:%d] cannot be processed.\n"
  230.        "Only [0:%d] vectors are valid",
  231.        x_re.q_lwb(), x_re.q_upb(), N-1);
  232.  
  233.   register int k;
  234.   register Complex * ak;
  235.  
  236.   for(ak=A,k=0; k < N; k +=8)
  237.   {
  238.     const double cu = W[ N/8 ].real();            // cos( pi/4 )
  239.  
  240.     register int i;
  241.     i = index_conversion_table[k];    // Index being "reverse" to k
  242.     double x0 = x_re(i);
  243.  
  244.     i = index_conversion_table[k+1];
  245.     double x1 = x_re(i);
  246.  
  247.     i = index_conversion_table[k+2];
  248.     double x2 = x_re(i);
  249.  
  250.     i = index_conversion_table[k+3];
  251.     double x3 = x_re(i);
  252.  
  253.     i = index_conversion_table[k+4];
  254.     double x4 = x_re(i);
  255.  
  256.     i = index_conversion_table[k+5];
  257.     double x5 = x_re(i);
  258.  
  259.     i = index_conversion_table[k+6];
  260.     double x6 = x_re(i);
  261.  
  262.     i = index_conversion_table[k+7];
  263.     double x7 = x_re(i);
  264.  
  265.     double t1 = x0 + x1;
  266.     double t2 = x2 + x3;
  267.     double t3 = x4 - x5;
  268.     double t4 = x4 + x5;
  269.     double t5 = x6 + x7;
  270.     double t6 = x6 - x7;
  271.     double t7 = t1 + t2;
  272.     double t8 = t4 + t5;
  273.  
  274. #define m2 t1
  275.     m2 -= t2;                // Forget t1 from now on
  276. #define m3 x0
  277.     m3 -= x1;                // Forget x0 from now on
  278.     double m7i = -cu*(t3 + t6);
  279. #define m4 t3                // Forget t3 from now on
  280.     m4 -= t6;
  281.     m4 *= cu;
  282.     double m5i = t5 - t4;
  283.     double m6i = x3 - x2;
  284.  
  285. #define s1 m3
  286.     double s2 = m3 - m4;        // Forget m3 from now on
  287.     s1 += m4;
  288. #define s3i m6i                // Forget m6 from now on
  289.     double s4i = m6i - m7i;
  290.     s3i += m7i;
  291.  
  292.     *ak++ = t7 + t8;                 // y0  
  293.     *ak++ = Complex(s1,s3i);            // y1
  294.     *ak++ = Complex(m2,m5i);            // y2
  295.     *ak++ = Complex(s2,-s4i);            // y3
  296.     *ak++ = t7 - t8;                // y4
  297.     *ak++ = Complex(s2,s4i);            // y5
  298.     *ak++ = Complex(m2,-m5i);            // y6
  299.     *ak++ = Complex(s1,-s3i);            // y7
  300.   }
  301.   assert( ak == A_end );
  302.   complete_transform();
  303. }
  304. #undef m2
  305. #undef m3
  306. #undef m4
  307. #undef s1
  308. #undef s3i
  309.  
  310. /*
  311.  *------------------------------------------------------------------------
  312.  *            Complex input with zero padding
  313.  *
  314.  * Note, if index j > N/2, its "inverse", index i is odd. Since the second
  315.  * half of the input data is zero (and isn't specified), x1, x3, x5, and x7
  316.  * in the Nussbaumer formulas above are zeros, and the formulas can be
  317.  * simplified as follows
  318.  *
  319.  *    t1 = x0                                   s1 = x0 + m4    y0 = t7 + t8
  320.  *    t2 = x2                                   s2 = x0 - m4    y1 = s1 + s3
  321.  *    t3 = x4            m2 = x0 - x2           s3 = m6 + m7    y2 = m2 + m5
  322.  *    t4 = x4            m3 = x0                s4 = m6 - m7    y3 = s2 - s4
  323.  *    t5 = x6            m4 = cu*(x4 - x6)                      y4 = t7 - t8
  324.  *    t6 = x6            m5 = -I * (x4 - x6)                    y5 = s2 + s4
  325.  *    t7 = x0 + x2       m6 = -I * x2                           y6 = m2 - m5
  326.  *    t8 = x4 + x6       m7 = -Isu* (x4 + x6)                   y7 = s1 - s3
  327.  *
  328.  *    cu = su = sin( pi/4 )
  329.  */
  330.  
  331. void FFT::input_pad0(
  332.          const Vector& x_re,    // [0:N/2-1] vector - Re part of input
  333.          const Vector& x_im)    // [0:N/2-1] vector - Im part of input
  334. {
  335.   are_compatible(x_re,x_im);
  336.   if( x_re.q_lwb() != 0 || x_re.q_upb() != N/2-1 )
  337.     _error("Sorry, vector [%d:%d] cannot be processed.\n"
  338.        "Zero padding is assumed, only [0:%d] vectors are valid",
  339.        x_re.q_lwb(), x_re.q_upb(), N/2-1);
  340.  
  341.   register int k;
  342.   register Complex * ak;
  343.  
  344.   for(ak=A,k=0; k < N; k +=8)
  345.   {
  346.     const double cu = W[ N/8 ].real();            // cos( pi/4 )
  347.  
  348.     register int i;
  349.     i = index_conversion_table[k];    // Index being "reverse" to k
  350.     Complex x0(x_re(i),x_im(i));
  351.  
  352.     i = index_conversion_table[k+2];
  353.     Complex x2(x_re(i),x_im(i));
  354.  
  355.     i = index_conversion_table[k+4];
  356.     Complex x4(x_re(i),x_im(i));
  357.  
  358.     i = index_conversion_table[k+6];
  359.     Complex x6(x_re(i),x_im(i));
  360.  
  361.     Complex t7 = x0 + x2;
  362.     Complex t9 = x4 - x6;
  363. #define t8 x4                // Forget x4 from now on
  364.     t8 += x6;
  365.  
  366.     Complex m2 = x0 - x2;
  367.     Complex m5(t9.imag(),-t9.real());
  368.     Complex m6(x2.imag(),-x2.real());
  369.     Complex m7(cu*t8.imag(),-cu*t8.real());
  370. #define m4 t9                // Forget t9 from now on
  371.     m4 *= cu;
  372.  
  373. #define s1 x0
  374.     Complex s2 = x0 - m4;        // Forget x0 from now on
  375.     s1 += m4;
  376. #define s3 m6                // Forget m6 from now on
  377.     Complex s4 = m6 - m7;
  378.     s3 += m7;
  379.  
  380.     *ak++ = t7 + t8;                 // y0  
  381.     *ak++ = s1 + s3;                // y1
  382.     *ak++ = m2 + m5;                // y2
  383.     *ak++ = s2 - s4;                // y3
  384.     *ak++ = t7 - t8;                // y4
  385.     *ak++ = s2 + s4;                // y5
  386.     *ak++ = m2 - m5;                // y6
  387.     *ak++ = s1 - s3;                // y7
  388.   }
  389.   assert( ak == A_end );
  390.   complete_transform();
  391. }
  392. #undef t8
  393. #undef m4
  394. #undef s1
  395. #undef s3
  396.  
  397. /*
  398.  *------------------------------------------------------------------------
  399.  *           Real input sequence with zero padding
  400.  * Again, since the input is a real sequence, some intermediate results
  401.  * are also real, that makes things simpler.
  402.  */
  403.  
  404. void FFT::input_pad0(
  405.          const Vector& x_re)    // [0:N/2-1] vector - Re part of input
  406. {
  407.   x_re.is_valid();
  408.   if( x_re.q_lwb() != 0 || x_re.q_upb() != N/2-1 )
  409.     _error("Sorry, vector [%d:%d] cannot be processed.\n"
  410.        "Zero padding is assumed, only [0:%d] vectors are valid",
  411.        x_re.q_lwb(), x_re.q_upb(), N/2-1);
  412.  
  413.   register int k;
  414.   register Complex * ak;
  415.  
  416.   for(ak=A,k=0; k < N; k +=8)
  417.   {
  418.     const double cu = W[ N/8 ].real();            // cos( pi/4 )
  419.  
  420.     register int i;
  421.     i = index_conversion_table[k];    // Index being "reverse" to k
  422.     double x0 = x_re(i);
  423.  
  424.     i = index_conversion_table[k+2];
  425.     double x2 = x_re(i);
  426.  
  427.     i = index_conversion_table[k+4];
  428.     double x4 = x_re(i);
  429.  
  430.     i = index_conversion_table[k+6];
  431.     double x6 = x_re(i);
  432.  
  433.     double t7 = x0 + x2;
  434.     double t9 = x4 - x6;
  435. #define t8 x4                // Forget x4 from now on
  436.     t8 += x6;
  437.  
  438.     double m2 = x0 - x2;
  439.     double m5i = -t9;
  440.     double m7i = -cu*t8;
  441. #define m4 t9                // Forget t9 from now on
  442.     m4 *= cu;
  443.  
  444. #define s1 x0
  445.     double s2 = x0 - m4;        // Forget x0 from now on
  446.     s1 += m4;
  447.     double s3i = -x2 + m7i;
  448.     double s4i = -x2 - m7i;
  449.  
  450.     *ak++ = t7 + t8;                 // y0  
  451.     *ak++ = Complex(s1,s3i);            // y1
  452.     *ak++ = Complex(m2,m5i);            // y2
  453.     *ak++ = Complex(s2,-s4i);            // y3
  454.     *ak++ = t7 - t8;                // y4
  455.     *ak++ = Complex(s2,s4i);            // y5
  456.     *ak++ = Complex(m2,-m5i);            // y6
  457.     *ak++ = Complex(s1,-s3i);            // y7
  458.   }
  459.   assert( ak == A_end );
  460.   complete_transform();
  461. }
  462. #undef t8
  463. #undef m4
  464. #undef s1
  465.